function [T ut ux P1 P2] = h2icu0(T, ux, uy, params)
  A = params.A;

  % given T ux, uy.

  % compute ut
  ut = sqrt(1+ux.^2+uy.^2);

  % compute spatial derivatives of initial conditions
  utdx = dx(ut);
  utdy = dy(ut);
  utdxdx = dx(utdx);
  utdxdy = dx(utdy);
  utdydy = dy(utdy);

  uxdx = dx(ux);
  uxdy = dy(ux);
  uxdxdx = dx(uxdx);
  uxdxdy = dx(uxdy);
  uxdydy = dy(uxdy);

  uydx = dx(uy);
  uydy = dy(uy);
  uydxdx = dx(uydx);
  uydxdy = dx(uydy);
  uydydy = dy(uydy);

  Tdx = dx(T);
  Tdy = dy(T);
  Tdxdx = dx(Tdx);
  Tdxdy = dx(Tdy);
  Tdydy = dy(Tdy);

  Tdt = 0;
  uxdt = (-1).*T.^(-1).*Tdx;
  uydt = (-1).*T.^(-1).*Tdy;
  utdt = ut.^(-1).*(ux.*uxdt+uy.*uydt);

  Tdtdx = dx(Tdt);
  Tdtdy = dy(Tdt);
  utdtdx = dx(utdt);
  utdtdy = dy(utdt);
  uxdtdx = dx(uxdt);
  uxdtdy = dy(uxdt);
  uydtdx = dx(uydt);
  uydtdy = dy(uydt);

  Tdtdt = (1/2).*(Tdxdx+Tdydy);
  uxdtdt = 0;
  uydtdt = 0;
  utdtdt = 0;

  P1 = (1/4).*(1+A).*T.^(-1).*(2.*Tdx.^2+(-2).*Tdy.^2+T.*((-1).*Tdxdx+Tdydy));
  P2 = (-1/2).*(1+A).*T.^(-1).*(T.*Tdxdy+(-2).*Tdx.*Tdy); 

  P1dx = dx(P1);
  P1dy = dy(P1);
  P2dx = dx(P2);
  P2dy = dy(P2);




end
